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1.  INTRODUCTION 

In  a recent  paper  [5],  Glassey  proposed  a special 
the  computation  o£  equilibrium  prices  in  a single-conmodlty  spatial  economic 
model.  The  model  was  formulated  as  a convex  quadratic  program  with  a 
separable  objective  function  and  with  certain  transshipewnt-type  constraints. 
Spatial  equilibrium  models  of  this  type  have  been  discussed  intensively  in 
the  literature . (See  Samuelson  [IS],  Takayama  and  Judge  [16]  and  more 
recently,  Folito  [14].)  Unfortunately,  the  original  algorithm  proposed  by 
Glassey  falls  to  work  in  certain  cases.* 

The  essential  purpose  of  this  paper  is  to  present  an  efficient  special- 
purpose  algorithm  for  solving  the  quadratic  program  arising  from  the 
aforementioned  economic  equilibrium  model.  The  proposed  algorithm  is  based 
on  an  efficient  implementation  procedure  of  the  parametric  principal  pivoting 
algorithm  for  solving  a parametric  linear  complementarity  problem.  (See  Cottle 
[2]  for  the  description  of  the  original  version  of  the  latter  algorithm  and 
Fang  [9]  for  that  of  its  implementation  procedure.)  This  latter  algorithm 
has  been  found  very  useful  in  various  other  contexts,  see  Kaneko  [9,10], 

Pang  [11,12]  and  Pang,  Kaneko  and  Hallman  [13]. 

The  reformulation  of  the  quadratic  program  as  a linear  complementarity 
problem  to  be  solved  by  the  parametric  principal  pivoting  algorithm  leads 
to  the  study  of  an  interesting  matrix  associated  with  a simple  digraph 
having  weights  on  the  nodes.  Several  basic  properties  of  such  a matrix  are 
derived.  These  properties  are  extremely  useful  in  the  development  of  the 
specialized  algorithm  for  solving  the  spatial  equilibrium  model. 

The  organization  of  the  paper  is  as  follows.  The  next  section  consists 

* In  a private  communication  [ 6 ] , Glassey  informed  us  that  the  statement 
of  the  algorithm  was  Incomplete. 


-purpose  algorithm  for 
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of  two  parts.  In  the  first  part,  we  review  the  formulation  of  the  spatial 
equilibrium  model  addressed  by  Glassey  and  derive  from  this  formulation  an 
equivalent  linear  complementarity  problem  with  a very  special  matrix  which 
is  intimately  related  to  the  underlying  network  structure  of  the  model. 

In  the  second  part,  we  review  the  parametric  principal  pivoting  algorithm 
and  its  implementation  procedure.  The  third  section  presents  some  basic 
properties  of  the  special  matrix  appearing  in  the  linear  complementarity 
formulation  of  the  economic  model.  In  the  fourth  section,  we  specialize 
the  parametric  principal  pivoting  algorithm  to  solve  this  linear  complementarity 
problem  and  show  how  the  algorithm  can  be  greatly  simplified  by  exploiting 
the  structure  of  the  matrix.  In  the  fifth  section,  ve  report  some  computational 
results  with  the  application  of  the  specialized  algorithm  in  solving  some 
randomly  generated  problems  of  considerably  large  size  and  draw  some  concluding 
remarks.  Finally,  in  an  appendix,  we  point  out  how  Glassey's  algorithm  can 
sometimes  fail  to  find  the  equilibrium  prices  and  present  a counterexample. 


I 


2.1  Model  Formulation.  The  problem  addressed  by  Glassey  [5]  is  to  find 


a certain  spatial  model 


Here  N denotes  the  total  number  of  regions  under  consideration,  p denotes 

9 

the  desired  equilibrium  price  in  the  or-th  region,  Ja  is  the  net  import  in 
that  region,  x is  the  export  from  region  a to  region  8 and  is  referred  to 


as  a flow  variable,  a is  the  (given)  equilibrium  price  in  the  absence  of 


imports  and  exports,  b is  a given  positive  number  which  is  related  to  the 


of  supply  and  demand,  finally,  c^  is  the  nonnegative  unit 
transportation  cost  from  region  a to  region  8 satisfying  the  triangle  inequality 


Incidentally,  the  above  model  is  a simplification  of  those 
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discussed  in  [16].  Due  to  its  simplicity,  efficient  special-purpose  algorithm 
can  be  developed.  As  noted  by  Glassey,  problem  (1)  forms  the  Kuhn-Tucker 
conditions  for  the  (convex)  quadratic  program 


minimize 


N , N 

^ + cag*«g  ' Va>  -ubJ*ct  ^ Ob)  and  (Id). 


Notice  that  by  (1b),  we  may  assume,  with  no  loss  of  generality,  that 
xaa  ■ 0 *•  Using  (la)  and  (1b)  to  eliminate  the  variables  p^  and  y^, 

we  may  easily  cast  problem  (1)  as  the  linear  complementarity  problem:  find 

z€rn(n-i)  9uch  ch4C 


q + Mz  > 0,  z > 0 and  zA(q  + Hz)  » 0. 


(2a) 


Here  q » (qfc)  and  z « (z^)  with  k denoting  the  ordered  pair  (a,0)  of  distinct 
indices  a and  0 are  N(N-1 ) -vectors  whose  k -components  are  defined  as 


"lc  + c«s 


and  z^  ■ x 


«3  * 


(2b) 


Moreover,  M ■ (m^)  with  k and  -C  denoting  certain  ordered  pairs  of  distinct 
indices  is  the  N(N-1)  x N(N-1)  matrix  whose  (k,£)-entry  is  defined  as 


b + b 
<»  9 


-<b*  + b9> 


mlU 


-b_ 


-b 


if  k ■ l ■ (a,0) 

if  k - (a, 0)  and  l • (3, a) 

if  k - (a,0),  l » (y,0)  end  a j*  y 

if  k - (a,3) . I - (ar.Y)  *nd  0 f*  y 

if  k - (a,0) , l - (0  ,y)  end  a + y 

if  k - (a,0),  l - (v,a)  and  0 d y 

otherwise. 


(2c) 
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The  matrix  M is  an  example  of  an  arc-arc  weighted  adjacency  matrix. 
Formally,  given  a simple  digraph  G - (V,A),  with  V and  A denoting  the  sets 
of  nodes  and  (directed)  arcs  of  the  graph  respectively,  and  a set  of  positive 
scalars  {b^g  representing  weights  on  the  nodes,  the  arc-arc  weighted 

adjacency  matrix  of  this  weighted  graph  G is  the  real  square  matrix  M such 
that,  associated  with  each  pair  of  arcs  k and  K,  in  A,  is  the  entry  defined 

as  in  (2c).  We  point  out  that  the  digraph  arising  from  the  spatial  economic 
equilibrium  model  is  complete. 

At  this  point,  we  should  emphasize  that  the  entire  development  of  this 
paper , including  all  the  results  and  the  proposed  algorithm  for  the  equilibrium 
model,  does  not  depend  on  the  completeness  of  the  simple  digraph  arising  from 
the  model.  Consequently,  our  approach  is  applicable  not  only  to  the  one 
treated  by  Glassey,  but  as  well  as  to  its  generalization  where  the  underlying 
digraph  is  arbitrary  and  is  not  necessarily  complete.  Nevertheless,  in 
the  remainder  of  this  paper  we  continue  to  assume  that  the  digraph  arising 
from  the  equilibrium  model  is  complete. 

The  proposition  below  identifies  the  first  property  of  an  arc-arc 
weighted  adjacency  matrix. 

Proposition  1 . Let  G ■ (V,A)  be  a simple  (weighted)  digraph  with  M denoting 
its  arc-arc  weighted  adjacency  matrix.  Then  M is  synmetric  and  positive 
sesd-def inite . 

Proof.  The  symmetry  is  obvious.  To  show  that  M is  positive  semi -definite, 
let  z - (zk)k g ^ be  an  arbitrary  vector.  By  an  easy  calculation,  we  may  deduce 

tSta.  - £ 

or€V 


where 


6 

y“  % fpa  ’ 0 6Sa  Z(a'0) 

with  P and  S denoting  the  sets  of  predecessors  and  successors  of  a 
Of  Of 

respectively.  Q.E.D. 

The  proof  of  the  proposition  suggests  that  an  arc-arc  weighted  adjacency 
matrix  may  be  related  to  the  node-arc  incidence  matrix  of  the  graph.  This 
relationship  is  made  explicit  in  the  next  proposition  whose  proof  is  easy  and 
thus  omitted. 

Proposition  2.  Let  G and  M be  as  in  proposition  1 . If  A is  the  node-arc 
incidence  matrix  of  the  digraph  G,  then 

M - ATDA 

where  D is  the  diagonal  matrix  whose  diagonal  entries  are  the  weights  b 

a 

of  the  nodes. 

As  we  shall  demonstrate  later,  M is  not  always  nonsingular,  or 
equivalently,  positive  definite. 

2.2.  Parametric  principal  Pivoting.  As  indicated  in  the  introduction, 
we  plan  to  solve  the  linear  complementarity  problem  (2)  by  the  parametric 
principal  pivoting  algorithm.  For  this  purpose,  we  briefly  review  this 
algorithm  and  its  implementation  procedure. 

Formally,  for  given  n-vectors  p and  q with  p > 0 and  nxn  matrix  M, 
the  par”1* trie  linear  complementarity  problem  is  to  find,  for  each  value 
of  the  parameter  \ lying  in  some  interval  [!,“),  an  n-vector  z satisfying 

q + Xp  + Mz  > 0,  z > 0 and  zT(q  + Xp  + Hz)  ■ 0.  (3) 
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Algorithm  1.  Parametric  Principal  Pivoting  Algorithm. 

Step  0.  (Initialization)  Let  the  initial  canonical  system  be  given  by 

w - q + Xp  + Mz 

where  q + Xp  > 0 for  sufficiently  large  values  of  X.  Put  j •*  0,  XQ  * • . 

Step  1 . (Ratio  test)  If  p < 0 (this  cannot  occur  initially) , let 
(w(X),z(X))  • (q  + Xp,  0)  for  X lying  in  the  interval  [X,Xj]  and  terminate. 
Otherwise,  determine  the  critical  index  r by 

xj+1  " ‘ Vpr  * max  f "Vpi  : > °1- 

Put  (v(X),z(X))  - (q  +Xp,0)  on  [X,  , ,X,  J. 

j+i  j 

— eP-  -•  (1x1  Diagonal  pivot)  If  m^  -0,  go  to  Step  3.  Otherwise, 

pivot  on  m^  and  let  w,z,p,q  and  M correspond  to  the  resulting  system. 

Replace  j with  j+1  aud  go  to  Step  1. 

Ste2_3.  (2x2  Block  pivot)  If  mlr  > 0 for  all  i,  redefine  X by  X - Xj+1 

and  terminate.  The  problem  (3)  has  no  solution  for  X < X.  Otherwise,  define 
the  critical  index  a by 


Perform  a pivot  operation  successively  on  m^  and  m#r  . Let  w,z,q,p  and  M 
correspond  to  the  resulting  system.  Replace  j with  j+1  and  go  to  Step  1. 

We  refer  to  Cottle  [1]  and  Graves  [71  for  a detailed  discussion  of  the 
theory  of  pivotal  algebra.  Under  the  condition  that  the  initial  M is 
a P-matrix  (i.e.  has  positive  principal  minors)  or  is  positive  semi -definite, 
the  algorithm  above  will  compute  a solution  function  z(X)  to  the  parametric 
linear  complementarity  problem  (3)  in  a finite  number  of  steps,  provided  that 
some  appropriate  tie-breaking  rule,  like  the  one  described  in  [7],  has  been 


s 


i 


[ 


t I 


e.  • 

' !* 


I 


incorporated  in  the  algorithm.  (See  [2].) 

The  2x2  block  pivots  described  in  Step  3 of  the  Algorithm  are  designed 
specially  for  positive  semi -definite  matrices.  They  are  redundant  if 
M is  a P-matrix.  In  [11],  the  first  author  has  proposed  an  efficient 
implementation  procedure  for  actually  carrying  out  the  1 x 1 diagonal  pivots. 
The  procedure  la  rephrased  in  the  algorithm  below. 

Algorlthm_2.  The  Parametric  Principal  Pivoting  Algorithm  Using  only 
Diagonal  pivots. 

Ste2_0.  Let  XQld  - Let  J - # and  I - {l,...,n). 

Ste£_1_.  Solve*  the  system  of  linear  equations  for  (q  , p ) : 

J «J 

MJJ(qJ,pP  " 


and  compute 


(vpi}  - <vv  - “iA’Pj5- 

Ste2_2.  If  P < 0 and  p > 0,  set  (z(X))_  - 0 and  (z(X))  « -q  - \p 
* J i J J J 

for  all  \ and  terminate.  Otherwise  determine  the  new  critical  value 

Xnew  * “*  f f - V’i  i 61},  max  { - : p < 0,  j € j} } 

and  let  k be  a maximizing  index.  Put  (z(X))_  - 0 and  (z(X))  - - q - Xp 
for  all  X in  the  interval  [Xnew,  XoldJ.  If  Xnaw  < X,  terminate.  Otherwise 
Mt  Xold  * Slew* 

* If  A is  an  nxm  matrix  and  I and  J are  subsets  of  {l,...,n}  and  {l,...,m} 
respectively,  then  by  A^j  we  mean  the  submatrix  of  A whose  rows  and  columns 

are  Indexed  by  I and  J respectively.  Similarly,  if  q is  an  n-vector,  then 
by  q^,  we  mean  the  subvector  of  q whose  components  are  Indexed  by  I. 


1 


! 


- 


■o— 
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Step  3.  Update  the  index  sets  I and  J: 


new 


new 


'JoldU  tkl 

££  k * Jold 

Jold  \ t k 1 

££  k 6 ■'old 

» • • • »nj\J  • 

new 

Go  to  Step  1. 

We  point  out  three  remarks.  First,  the  parametric  principal  pivoting 
algorithm  can  be  easily  used  to  solve  a linear  complementarity  problem 


of  the  form  (2a).  In  fact,  it  suffices  to  choose  1*0.  Second,  Algorithm 
2 is  also  applicable  when  M is  positive  semi -definite,  provided 
that  it  is  not  necessary  to  perform  the  2x2  block  pivots.  As  a matter 
of  fact,  we  shall  show  in  Section  4 that  this  is  indeed  the  case  when  the 
linear  complementarity  problem  (2)  is  solved  by  Algorithm  1.  Finally,  if 
a linear  complementarity  problem  with  a positive  semi -definite  matrix 
is  solved  by  Algorithm  1 and  if  all  the  pivots  are  1x1  diagonal,  then 
the  algorithm  must  terminate  with  a solution  to  the  problem. 


<3 


^ • BASIC  PROPERTIES  OF  AN  ARC -ARC  WEIGHTED  ADJACENCY  MATPTY 

In  this  section,  we  establish  some  basic  properties  of  an  arc-arc 
weighted  adjacency  matrix  M associated  with  a general  simple  digraph  having 
weights  on  its  nodes.  Our  objectives  are,  first,  to  characterize  the 
nonsingularity  of  an  arbitrary  principal  submatrix  of  M in  terms  of  the 
structure  of  the  "associated  subgraph"  of  the  submatrix,  and  then,  to  show 
how  the  solution  of  the  various  systems  of  linear  equations  required  in 
Step  1 of  Algorithm  2 can  be  achieved  by  identifying  the  decendents 
of  each  node  in  the  associated  subgraph  of  MTT. 

Noticing  that  each  principal  submatrix  of  an  arc-arc  weighted  adjacency 
matrix  is  itself  an  arc-arc  weighted  adjacency  matrix  associated  with  an 
obvious  subgraph,  we  shall  develop  the  results  in  terms  of  an  arc -arc  weighted 
adjacency  matrix  M defined  with  respect  to  a fixed  but  arbitrary  weighted 
digraph  G with  the  understanding  that  these  results  apply  readily  to  each 
of  the  principal  submatrices  of  possibly  another  arc-arc  weighted  adjacency 
matrix.  In  particular,  this  digraph  G is  not  necessarily  the  one  arising 
from  the  spatial  equilibrium  model  discussed  in  the  last  section. 

Recall  that  a cycle  in  a digraph  is  a sequence  of  arcs  connecting  a 
node  to  itself.  The  directions  of  the  arcs  are  irrelevant  in  the  cycle.  A 
minimal  cycle  is  one  which  contains  the  minimum  number  of  arcs.  It  is 
important  to  observe  that  if  an  arc  is  incident  to  two  nodes  in  a minimal 
cycle,  then  the  arc  must  be  one  of  those  in  the  cycle.  Throughout  the 
paper,  all  digraphs  are  simple. 


Figure  2:  The  outer  cycle  is  not  minimal 


Propoaltion  3.  If  the  digraph  G contains  a cycle,  then  its  associated 
arc-arc  weighted  adjacency  matrix  M is  singular. 

Proof.  It  suffices  to  exhibit  a nonzero  vector  v such  that  Mv  ■ 0.  Choose 
a minimal  cycle  in  G and  let  A1 , , . . . ,A^  be  the  arcs  in  the  cycle.  Fix 


follows 


if  h • A^  for  some  i and  if  A^  and  A^  are  oriented 
in  the  same  direction 
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We  claim  chac  Mv  « 0.  Consider  a specific  component,  say  (Mv)^  with  4 denoting 

a certain  arc.  If  the  arc  4 is  not  connected  to  any  node  in  the  minimal  cycle, 
then  obviously  (Mv)^  - 0.  On  the  other  hand,  if  4 is  incident  to  one  single 

node  in  the  cycle,  say  node  a,  then  the  4-row  of  M is  given  by 

*1  A2  ...  At  ...  At,  ...  Al-1 

[ ...  0 0 ...  ...  #2ba  •••  0 0 ...  ] 

where  A and  A , are  the  two  arcs  in  the  cycle  which  are  incident  to  node  a 

U U 


and  where  Sj(s2)  1*  equal  to  +1  depending  on  whether  the  arcs  4 and  A^A^, ) 

are  oriented  in  the  opposite  or  the  same  direction.  If  both  pairs  of  arcs 
(4  and  A£)  and  (4  and  At,)  are  oriented  in  the  same  way  (i.e.,  if  s^  ■ s^), 
then  the  arcs  Afc  and  A£l  must  be  oriented  in  the  opposite  direction  (i.e., 

then  v » - v );  on  the  other  hand,  if  the  two  pairs  of  arcs  (4  and  A ) 
\ t 


and  (4  and  A , ) are  oriented  differently  (i.e.,  if  s,  - - s0) , then  A and  A . 

c 1 4 c t 

must  be  oriented  in  the  same  direction  (i.e.,  then  v.  * v.  ).  Consequently, 

t t' 


in  either  case , we  have 


(M»)t  - 


0 . 


Finally,  if  4 is  incident  to  two  nodes  in  the  cycle,  then  it  must  coincide 
with  one  of  the  L arcs  by  minimality  of  the  cycle.  Let's  suppose 
4 ■ A ■ (a, 9).  Then  the  4-row  of  M is  given  by 


t 


A2  ...  Afc  ...  A ...  A ...  A^ 
0 ...(ba  + b0)  ...  s1btf  ...s2b0  ...  0 


0 ...  ] 


where  A^  and  A^  are  the  two  arcs  adjacent  to  At  In  the  cycle  and  where  s^sp 

is  equal  to  ±1  depending  on  whether  A and  A (A  ) are  oriented  in  the 

t u,  u2 

opposite  or  the  same  direction.  By  using  a similar  argument,  we  may  deduce 


' <b«  + Vv  + *ib«\  + *2Va  ' 0 • 

U1  u2 


Consequently,  M is  in  fact  singular.  Q.E.D. 


Given  a nonsingular  principal  submatrix  of  a square  matrix  M,  there 
exists  a permutation  matrix  P such  that 


T 

P HP 


M11 

M12 

*21 

«22 

The  Schur  complement  of  in  M is  the  matrix  M22  - Moreover, 

the  Schur  determinental  formula  (see  Cottle  [3]  e.g.)  is 


det  M - det  Mn  x det  (M^  - 


The  next  proposition  is  a partial  converse  of  the  last  one. 

Proposition  4.  if  the  digraph  G is  a tree,  then  the  arc-arc  weighted 
adjacency  matrix  M is  nonsingular. 

Proof.  We  use  induction  on  n(>  2)  the  number  of  nodes  in  G.  The  assertion 
is  obvious  for  n - 2.  Suppose  that  it  is  true  for  a tree  with  n nodes. 
Consider  a tree  G with  n + 1 nodes.  Since  G is  a tree,  it  contains  an  end 
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node,  say  a . Let  A be  the  only  arc  in  G incident  to  a and  let  9 be  the  other 

end  point  of  the  arc  A.  Let  G'  be  the  subtree  obtained  from  G by  deleting 

the  node  a and  the  arc  A.  Assign  the  same  set  of  weights  to  all  nodes 

in  G1  except  for  the  node  9 where  we  assign  the  weight  1/(b_1  + b"1).  Let 

at  p 

M'  be  the  arc-arc  weighted  adjacency  matrix  obtained  from  the  (weighted) 
digraph  G' . By  induction  hypothesis,  M'  is  nonsingular.  We  now  show  that 
M'  is  precisely  the  Schur  complement  of  the  diagonal  entry  MAA  in  the 

original  matrix  M.  With  no  loss  of  generality,  we  may  assume  that 


1 maa 

MA2 

M - 

' “2A 

**22 

where  is  given  by 

MA 2 " t0  0 b9  •**  b9  " b9  " b9  0 •“  0 1 

with  the  plus  (minus)  b^  appearing  in  those  arcs  (i.e.,  columns)  which 

are  incident  to  the  node  9 and  oriented  in  the  opposite  (same)  direction 
as  the  arc  A.  Moreover,  the  principal  submatrix  of  corresponding  to 

those  arcs  which  are  incident  to  the  node  9 is  given  by 


where  E is  the  matrix  of  1's  and  where  £ is  the  diagonal  matrix  whose 
diagonal  entries  art  the  weights  of  those  nodes  (except  for  or)  which  are 
adjacent  to  the  node  9 • The  same  principal  submatrix  of  -M91MAAK,g  is 
given  by 


• S r E 1 ~E] 

, + b,  L-e  1 iJ  . 


— - 
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la  face,  all  entries  of  the  matrix  -M21M^M12  are  zero  except  for  those  of 
this  principal  submatrix.  it  then  follows  that  all  entries  of  the  Schur 

complement  M22‘M21MAA.M12  are  the  aaae  aa  the  corresponding  ones  of  M^,  except 
for  those  in  the  principal  submatrix  which  is  equal  to 


~ E 

Lp 

i r 1 

~ E 

Lll 

[ 6 * + b»  ) 

-E 

b"1  + b"1 
a p 

-E 

It  18  ohrtou.  th«t  Mj2  . . M'  . Since  both  end  It. 

Schur  complement  are  nonsingular,  it  follows  from  the  Schur  determinencal 
formula  that  so  is  M.  This  completes  the  inductive  step  and  the  proof.  Q.E.D. 
Remark.  The  above  two  propositions  can  also  be  proved  by  using  Proposition  2 
and  the  tree-property  of  the  basis  matrix  of  a linear  transhipment  problem 
(see  Dantzig  [4]). 

Corollary  5.  If  the  digraph  G is  a forest,  then  the  arc-arc  weighted 
adjacency  matrix  M is  nonsingular. 

Proof:  It  is  sufficient  to  observe  chat  such  a matrix  is  block  diagonal 

with  each  diagonal  block  corresponding  to  a tree.  The  corollary  then 
follows  immediately  from  Proposition  4.  Q.  E.  D. 


Corollary  6.  Let  M be  the  arc-arc  weighted  adjacency  matrix  associated 
with  the  (weighted)  tree  G.  Then 

det  m - (»  > x (IlT1) 

where  both  the  product  and  the  sunmatlon  range  over  all  the  nodes  in  G. 

Proof.  We  use  induction  on  the  number  n of  nodes  in  G.  The  formula  is 
certainly  correct  for  n - 2.  Suppose  that  it  is  true  for  a tree  with  n nodes. 
Consider  now  a tree  G with  n+1  nodes.  Using  the  same  notations  as  in  the 
proof  of  Proposition  4,  we  have  by  the  Schur  determinencal  formula. 


det  M ■ (b  + b ) x det  M’ 

a p 


Now  che  indueclon  hypothesis  implies  chee 
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det  M'  - 0 b x I/O*"1  + b'1)  x £ b’1 


Y*»,9 


« 9 


ell  y 


Pram  these  two  equations,  the  desired  formula  tor  det  M follows  readily. 

Q.  E.  0. 

Combining  Corollary  5 sol  Proposition  3.  we  obtain  the  theorem  below 
which  characterises  the  turns  lngularlty  of  an  arc-arc  weighted  adjacency 
matrix  associated  with  a (weighted)  digraph. 

Theorem  7.  Let  G be  a digraph  with  positive  weights  on  the  nodes.  Then 
the  arc-arc  weighted  adjacency  matrix  is  nonsingular  if  and  only  if  the 
graph  G is  a forest. 

Let  A - (or, 3)  be  a fixed  but  arbitrary  arc  in  a tree  G.  We  say 
that  a node  y which  is  different  from  or  and  9,  is  on  level  1 with  respect  to 
the  arc  A if  y is  incident  to  A.  Inductively,  given  a node  y not  on  level 
k - 2,  we  say  chat  it  is  on  level  k with  respect  to  the  arc  A if  it  is 
adjacent  to  one  of  the  nodes  on  level  k - 1.  The  nodes  at  and  % are  considered 
to  be  on  level  0.  Given  two  nodes  y,  and  y£  on  levels  k^  and  k2  respectively 

with  k1  < k2,  we  say  that  y2  i«  « descendant  of  y,  If  the  nodes  In  the 
unique  path  connecting  node  y^  to  node  y2  are  on  strictly  descending  levels. 

By  convention,  we  consider  a node  as  a descendant  of  itself.  We 
shall  denote  by  the  set  of  descendants  of  the  node  y.  Notice  that 

Da  sad  D.  partition  the  set  of  nodes  in  G.  We  call  the  number  E b“* 

a €d  ff 


» i 
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modified  weight  of  the  node  y,  and  denote  it  by  b . 


© © 

/ v/v 

© bo  o 

/ \ 

3)  n 


level  0 


level  1 


level  2 


level  3 


Figure  3:  Different  levels  in  a tree 
Note:  The  node  y^  is  a decendent  of 
y^  but  not  of  y^  . 


The  two  propositions  below  are  fundamental  to  the  solution  of  a system 
of  linear  equations  involving  an  arc-arc  weighted  adjacency  matrix. 
Proposition  a.  Let  G be  a weighted  tree  with  M being  its  arc-arc  weighted 
adjacency  matrix.  Let  A - (er,9)  be  any  arc  in  G.  Let  be  the  A-column 
of  M with  the  diagonal  entry  deleted  and  let  M'  be  the  principal. submatrlx 

of  M with  both  the  A-row  and  the  A-column  deleted.  Then  the  vector  (M1 ) 

is  given  by  the  following:  if  k * A is  the  arc  (5,y)  with  v€D,,  then 

0 

fS  /sum  if  1c  and  A are  oriented  in  the  opposite  direction 


-b  /sum  if  k and  A are  oriented  in  the  same  direction 
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where 


if  y€D„ 


if  Y € 0^  . 


Proof.  The  matrix  M*  is  block  diagonal  with  the  two  diagonal  blocks  Mj 

and  being  the  principal  subnatrlces  of  M*  whose  arcs  are  incident  to  nodes 
in  D#  and  respectively.  Due  to  this  block  diagonal  structure  of  M',  we 

may,  with  no  loss  of  generality,  assume  that  the  arc  A is  incident  to  an 
end  node,  say  8 ; or  in  other  words,  we  assume  that  M£  is  vacuous.  Let  G' 

be  the  subtree  obtained  from  G by  deleting  the  arc  A and  the  node  3 . We 
may  further  assume  that  M'  has  the  form 


Mi1 

Mi2 

«21 

“22 

where  MJ  ^ is  the  principal  submatrix  of  M'  whose  arcs  are  incident  to  a. 
With  respect  to  this  partitioning,  we  may  write 


* b~ 
1A  a 


where  is  a vector  of  plus  and  minus  ones  such  that  (M^)^  ■ 1 (- 1 ) if  the 

arcs  i and  A are  oriented  In  the  opposite  (same)  direction.  If  we  write 
(M')  m/xA  according,  to  the  above  partitioning  of  M',  we  may  easily  deduce 


(Mj,  . 


19 


By  the  proof  of  Proposition  4 and  an  inductive  argument,  it  is  not  difficult 
to  show  that  the  Schur  complement  (MJ , - is  the  arc-arc 

weighted  adjacency  matrix  associated  with  the  subtree  G"  obtained  from  G’ 
by  deleting  all  the  arcs  not  Incident  to  the  node  a and  all  the  nodes  not 
adjacent  to  a.  The  weight  of  a node  a in  the  subtree  G"  Is  equal  to  (b^)” ' 

If  a + a and  equal  to  btf  if  a - cr.  In  fact,  this  Schur  complement  Is 
given  explicitly  by 

Mii  ' ' 1 + W»)' 

where  £ is  the  diagonal  matrix  whose  diagonal  entries  are  the  inverses  of  the 
modified  weights  of  the  nodes  adjacent  to  a in  G".  According  to  Sherman- 
Morr is on -Wood bury  formula  (see  Householder  [8,  p.  124]  e.g.),  we  obtain 


‘i  - b«  1 


ia 


which  gives 

*1  * E''MiA/S«  • 

This  establishes  the  desired  formula  for  ((M')-1M1A)k  provided  that  k is 

incident  to  node  a.  To  establish  the  formula  for  k not  incident  to  a,  we 
evaluate  It  is  not  difficult  to  see  that  M£2  is  again  a block  diagonal 
matrix  with  each  diagonal  block  being  the  principal  submatrix  whose  arcs 
are  incident  to  nodes  in  for  some  node  ff  i*  a in  the  tree  G".  Again  we  may 
assume  that  M£2  consists  of  just  one  single  block  with  all  the  arcs  incident 
to  a certain  descendant  of  a specific  note  oj*ar  in  G".  We  then  have 


where  A'  is  Che  ere  connecting  cr  end  a,  sey  A'  « (a, a),  end  M^,  is  defined 
1$  e feshien  similar  to  ebove.  By  decomposing  x2  mf  *21  \ *21 

» X-_/ 
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corresponding  Co  Chose  srcs  incident  Co  Che  node  a end  using  Che  seme 
argument,  ve  may  easily  deduce 

where  £'  is  Che  diagonal  matrix  whose  diagonal  entries  are  Che  inverses  of 
Che  modified  weights  of  Che  nodes  adjacent  to  a.  Substituting  xag  Into  Che 
above  expression  for  Xj^  immediately  yields  the  desired  formula  for 

for  k incident  to  o.  The  proof  of  the  proposition  can  now 
be  completed  by  repeatedly  using  the  above  argument  to  all  branches  of 
the  tree  G.  Q.  E.  D. 

Proposition  9.  Let  G,M,A,M^  and  M'  be  as  in  Proposition  8.  Then 

“aa  - <Xia>V>-'m,a  ■ (S/1  + <V" 

Proof.  To  simplify  the  notations,  we  assume,  as  we  did  earlier,  that 
the  arc  A is  incident  to  the  end  node  B . Then  by  Corollary  6 and  the 
Schur  deteminental  formula,  we  have 


w 


« det  M/det  M' 


V - <«,a>V>‘'m,a 


-1 


Y 

y*9 


VbY 


■ + <V'' 


«•  desired.  Q.  E.  D. 
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Using  the  same  notations  as  in  Proposition  8,  we  consider  the  solution 
of  the  system  of  linear  equations  Mx  - d for  some  given  vector  d.  We  may 
assume  with  no  loss  of  generality  that  M'  is  a leading  principal  submatrix 
of  M.  Partitioning  the  vectors  x and  d accordingly  we  may  write 


which  Implies 

(MAA  ’ (M1A)TcM,)’1m1A)xA  " dA  • (M1A)T(M')’,d*  (4a) 

and 


x'  - (M'rV  - (M,)”1M1A*A  . 


(4b) 


Combined  with  Propositions  8 and  9,  these  latter  two  formulas  show  how  the 
two  vectors  x and  (M1)  'd'  can  be  efficiently  computed  from  each  other.  In 
fact,  suppose  that  (M')-,d'  is  known,  then  Proposition  8 shows  that  (M1A)T(M')_1 
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can  be  obtained  by  identifying  the  descendants  of  each  node.  Together 

with  Proposition  9,  formula  (4a)  yields  the  component  x^  readily.  Substituting 

XA  and  using  Proposition  8 again,  we  may  easily  compute  x'  from  (4b).  Conversely, 

if  x is  known,  then  (M*)  *d'  can  be  obtained  readily  from  (4b)  as  well. 
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4.  SPECIALIZATION  OF  THE  PARAMETRIC  PRINCIPAL  PIVOTING  ALGORITHM 


In  this  section,  ve  specialize  the  parametric  principal  pivoting 
algorithm  to  solve  the  linear  complementarity  problem  (2)  arising  from  the 
spatial  equilibrium  model.  Our  purposes  are  (1)  to  show  that  the  problem 
can  be  solved  by  perfoxmlng  the  1 z 1 diagonal  pivots  exclusively  and  (2)  to 
derive  from  this  result  and  those  established  in  the  last  section  an  efficient 
specialized  algorithm  for  solving  the  problem.  To  achieve  these,  ve  first 
state  two  preliminary  results  having  to  do  with  the  application  of  the 
parametric  principal  pivoting  algorithm  for  solving  a general  parametric 

| 

linear  complementarity  problem  with  a symmetric  positive  semi -definite  matrix. 

Lemma  10.  Let  M be  any  symmetric  positive  semi -definite  matrix.  Consider 
the  solution  of  the  parametric  linear  complementarity  problem  (3)  by 
Algorithm  1 . Let  J denote  the  index  set  of  the  basic  z-varlables  at  each 
iteration.  Then  Hjj  is  nonsingular.  Moreover,  if  a 2x2  block  pivot  occurs 
at  a certain  iteration  and  if  r and  s are  the  two  critical  indices  obtained  in 
Steps  1 and  3 of  the  algorithm,  then  rf  J and  s 6 J. 

The  two  assertions  contained  in  this  Lemma  are  rather  well-known  in 
the  theory  of  pivotal  algebra.  Their  proofs  are  thus  omitted.  (The 
authors  are  grateful  to  Professor  I.  Kaneko  for  pointing  out  this  fact.) 

Corollary  1 1 . Let  I be  the  index  set  of  the  currently  nonbasic  z-varlables. 

If  the  twag-mii—  ratio  in  Step  1 of  Algorithm  1 does  not  occur  at  the 
nonbasic  index  l€l,  then  the  variable  z^  can  not  become  basic  at  the  next 
pivot. 

Proof.  Under  the  given  assumption,  the  only  possible  way  for  zi  to  become 
basic  is  for  i equal  to  the  critical  index  s.  But  Lemma  10  asserts  that  this 
is  impossible.  Q.  E.  0. 
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l 
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With  these  two  general  results,  we  proceed  to  solve  the  linear 
complementarity  problem  (2).  Recall  that  the  lower  limit  X of  the  parameter  X 
is  zero.  Combining  Lemma  10  with  Theorem  7,  we  obtain 

Theorem  12.  Throughout  the  solution  process,  the  set  of  basic  flow  variables 
can  not  contain  a cycle. 

In  [4],  Glassey  showed  that  a solution  to  (1)  can  be  found  which 
contains  no  transshipment  of  flows . The  theorem  below  extends  this  result. 

Theorem  13 . Let  the  parametric  vector  p be  chosen  such  chat  all  components 
are  equal.  Then  throughout  the  solution  process,  there  is  no  transshipment 
of  flows. 

We  need  three  lemmas  to  establish  the  theorem.  From  now  on,  we  let 
the  vector  p be  chosen  as  stated  in  the  theorem. 

Lemma  14.  Let  J be  the  index  set  of  the  currently  basic  flow  variables. 

Suppose  that  k • (or, 3)  €J.  Then  the  nonbasic  variable  z^  with  L ■ (g.ar) 

\ 

can  not  become  basic  at  the  next  iteration. 

Proof.  Notice  that  Theorem  12  implies  that  L£J.  The  ^-component  of 
the  current  q-vector  before  the  next  pivot  is  given  by 

h ' \ - • 


It  is  easy  to  verify  that 

'-1 

1, 


Wj  - 


Therefore,  we  have 


if  j - k 
otherwise  . 


*1  m ’t  + \ ' c»3  + e».  2 0 ' 
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Hence,  according  to  Corollary  11  and  the  termination  criteria,  the 
variable  z^  can  not  become  baaic  at  the  next  iteration.  q.  e.  D. 


Lemma  15.  Let  J be  aa  in  Lemma  14.  Suppose  that  kj  - (or, 3)  and 
^2  " (®*Y)  with  p + y **•  ia  Then  :he  nonbasic  variables  * and  z 

h £ 

where  - (3,y)  and  l2  m (YiP)  can  not  become  basic  at  the  next  iteration. 
Proof.  Notice  that  Theorem  12  implies  that  i J and  -C^gj.  The 

-tj -component  of  the  current  q -vector  is  given  by 


\ - % - v(Mjj)  qj  • 


It  is  not  difficult  to  show  that  (cf.  the  proof  of  Proposition  3 e.g.) 


if  j - k. 


if  J - k. 


otherwise. 


Consequently,  we  have 

\ • \ + {\  • V ■ C°8+  cSv  • C«Y  2 0 

where  the  last  inequality  follows  from  the  triangle  inequality  of  the 
transportation  costs.  Hence  according  to  Corollary  11  and  the  termination 
criteria,  the  variable  z^  can  not  become  basic  at  the  next  iteration. 


Similarly,  we  may  establish  the  same  conclusion  for  z 


Q.  E.  D. 
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Lemma  16.  Let  J be  as  in  Lemma  14.  Suppose  that  J contains  no 
transshipment  of  flows.  If  the  arc  -L  is  such  that  JU  [t]  contains  a cycle, 
then  the  nonbasic  variable  z^  can  not  become  basic  at  the  next  iteration. 

In  fact,  the  L -component  of  the  current  q-vector  is  nonnegative. 

Proof.  Theorem  12  implies  that  J contains  no  cycle  and  that  -L  C J.  Suppose 
that  the  cycle  In  J U { t } contains  an  even  number  of  arcs . By  the  above 
two  lemmas,  we  may  assume  that  this  number  Is  at  least  four.  Notice  that  t 
must  be  one  of  these  arcs  and  that  there  are  precisely  two  arcs  adjacent 
to  each  arc  In  the  cycle.  Let  L » (0 ,y)  and  let  Che  two  arcs  adjacent 
to  i be  k1  ■ (Oj.fij)  and  k^  “ Becauae  the  even  number,  there 

are  two  cases:  (0  « and  y • ^ or  (Y  * «1  and  0 - 62>.  Consider  the 
first  case,  namely  0 » af^  and  y « 5^  (cf.  Figure  4).  The  ^-components  of 
the  current  q-  and  p -vector  are  given  by 


Figure  4.  An  even  cycle  with  no  transshipment 
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Vjj’j  “d  h ■ - vita  • 


By  Che  proof  of  Proposition  3 , one  can  easily  deduce  chat 


<vta  j 


f: 


If  j Is  not  In  Che  cycle 

If  j Is  In  Che  cycle  and  oriented  In  the 
opposite  direction  as  t 

If  j Is  In  the  cycle  and  oriented  in  the 
same  direction  as  1 . 


From  this  and  an  easy  calculation,  we  may  deduce 


q,  - E c -Sc 

L j j 1 1 


and  p.  ■ 0 

V 


where  the  first  and  second  sums  in  the  expression  for  q . range  over  those 
arcs  in  the  cycle  which  are  oriented  in  the  same  and  opposite  directions 
as  l respectively.  Notice  that  p^  - 0 follows  from  the  fact  that  the 
cycle  contains  an  even  number  of  arcs.  Consequently,  according  to  the 
ratio  test  in  Step  1 of  Algorithm  1 and  Corollary  11,  will  not  become 
basic  at  the  next  pivot.  Moreover,  we  must  have  q > 0 because  the  non- 

v “ 

negativity  of  the  component  q^  + Xp^  is  maintained  throughout  the  algorithm. 

Consider  the  second  case,  namely  y » and  0 - 82  (cf.  Figure  5).  By 
the  same  token,  we  may  similarly  deduce 


\ S cj  ‘2c 
J J i 1 


where  the  two  summations  are  ranging  over  the  same  seta  of  arcs  as  before. 
By  the  triangle  inequality,  it  follows  that 


1 


Figure  5.  An  even  cycle  with  transshipment 

^ * (®2 » ® j ) • Xk*  last  Inequality  follows  from  the 

^isst  case.  Consequently , the  variable  can  not  become  basic  at  the 
next  pivot.  Similarly,  we  may  establish  the  same  assertion  in  the  case  where 
the  cycle  is  odd.  Q.  g.  d. 

Proof  of  Theorem  73.  We  use  induction  on  the  number  of  iterations.  The 
assertion  is  obviously  true  at  the  first  iteration.  Suppose  that  it  is  true 
after  t iterations  and  the  algorithm  is  entering  its  (t  + 1)st  iteration. 

Let  J be  the  index  set  of  the  currently  basic  flow  variables.  By  induction 
hypothesis  J contains  no  transshipment  of  flows.  Moreover,  by  Theorem  t2,  J 
contains  no  cycle  as  well.  It  suffices  to  show  if  k - (er,(J)  € J,  then  all 
variables  z ^ with  i - (0,y)  or  t ■ (y,a)  can  not  become  basic  at  the  next 


pivot.  To  be  specific,  we  let  4 » (0  ,y) . By  Lemma  14,  we  may  assume  that 
Y + or.  Moreover,  by  Lemma  15,  we  may  assume  that  i - (ar,y)  £ J.  Finally 
by  Lemma  16,  we  may  assume  that  JU  (4}  contains  no  cycle  as  well.  The 
current  i-  and  4-components  of  the  p -vector  are  given  by 


Pt  • Pt  - “ijHjjPj  *nd  ft  " p4  " MtJMJJPJ 


Subtracting  the  two  equations  gives 


h ■ n • <Mu  - 


■ V1**  • \ 


Similarly,  we  may  deduce 


V<cSy  + c.9,s°- 


Hence,  it  follows  that 


- q^Pt  + - (-q4  + VPt  + qtpk 


Suppose  that  pi  < 0.  If  the  maximum  ratio  occurs  at  the  arc  4 such 
that  - q^  / P^  > 0 » then  we  have 

0 < - \ \ 


so  that  + Xpt  < 0 for  nonnegative  values  of  X.  But  this  contradicts 
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the  fact  that  q^  + Xp^  is  kept  nonnegative  throughout  the  algorithm. 
®iar®^ore»  ^ S 0 » then  either  the  algorithm  terminates  with  the 
arc  4 remaining  nonbasic  or  else  the  maximum  ratio  cannot  occur  at  the 
•rc  4.  In  either  case,  the  theorem  is  proved. 

Suppose  now  pt  > 0.  Then  we  have 

- Vi  + Vt  s Vk  • 

Now,  if  the  left-hand  term  is  nonnegative,  then  so  is  which  then  implies 
that  q^  > 0.  In  this  case  the  variable  z ^ can  not  become  basic 
in  the  next  iteration.  On  Che  other  hand,  if  the  left-hand  term  is 
negative,  then  we  would  have 

- <<4  '*0  < ’ («i  7»i> 

unless  p^  < 0.  Consequently,  the  maximum  ratio  will  not  occur  at  the 
arc  4.  The  desired  conclusion  therefore  follows  from  Corollary  11. 

Q.  E.  D. 

We  may  now  state  and  prove  our  principal  result. 

Theorem  17.  Let  the  parametric  vector  p be  chosen  to  be  the  vector  of 
ones.  Then  throughout  the  solution  process,  all  pivots  are  1x1  diagonal. 

In  other  words,  all  the  diagonal  pivot  entries  are  nonzero. 

Proof.  We  use  induction  on  the  number  of  pivot  iterations.  The  assertion 
is  certainly  true  at  the  first  pivot  because  all  diagonal  entries  of  M are 
positive.  Suppose  that  the  assertion  is  true  for  t iterations  and  the 
algorithm  is  entering  its  (t  + 1)st  iteration.  Let  J be  the  index  set  of 
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k**^c  variables  at  the  and  of  the  t-ch  iteration.  Suppose  now  that  the 
maxima  ratio  occurs  at  the  critical  index  r such  that  the  corresponding 
diagonal  pivot  entry  is  zero.  Then  rfj  and  the  matrix 


is  singular.  Consequently  the  set  of  arcs  in  JU{r}  must  contain  a cycle 
by  Theorem  7.  Since  J contains  no  transshipment  of  flows.  Lemma  16  implies 
that  the  current  critical  value  of  X is  nonpositive.  In  other  words,  the 
algorithm  will  terminate  without  performing  further  pivots.  Q.  E.  D. 


J 


: ■ 
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Sunmarlzing  the  discussion,  we  now  present  the  promised  algorithm 
for  solving  the  linear  complementarity  problem  (2)  arising  from  the 
spatial  equilibrium  model.  The  algorithm  is  a straightforward  specialization 
of  Algorithm  2 to  this  application. 


The  Algori£hmT 

Step  0 (Initialization)  Let  J - 0.  Determine 


X - max  ( - (ea  - ap  + caJ  ) : or  0 , « , 0 €V  ] 


If  X < 0,  terminate.  A solution  is  given  by  x - 0.  Otherwise,  let 
It  * (<*,0)  be  a maximizing  arc.  Let 


J’  - fj 


if  k £ J 
if  k 6 J . 


Step  1 (Search  for  descendants)  Identify  the  descendants  of  the  nodes 
or  and  0 which  are  incident  to  arcs  in  J'.  Determine  the  modified  weights 
of  these  descendants  (including  the  nodes  or  and  0 as  well).  Compute  the 


vector  Hj!  K according  to  Proposition  7*. 


Step  2 (Updating  of  basic  components)  Set 

-1  -1 


(*k  * V " [ (*k  * V / + 60  1 J 

(p^  > 5jj) 


if  k (£  J 
if  k 6 J 


* Die  set  J'  may  be  a forest  but  not  a tree.  The  proposition  is  nevertheless 


applicable  because  obviously,  (m"! j.MJflc)j  - 0 if  the  arc  J is  not  connected 

to  k.  Hence  we  need  to  compute  only  those  components  which  are  connected 
to  k. 


i 1 


J 


(Pji  » * (Pj*  * ” **J' j'^j'k^k  * ^k^ 


Step  3 (Updating  of  basic  index  set)  Set 


J-  fj  U { k ) 

1 

J\£k} 


if  k t J 


if  k € J 


Step  4 (Updating  of  nonbasic  components)  Lee  I • JC  . Compute 
(Pj  » 5j)  • (Pj  » *j)  - MjjCPj  » 5j)  • 


Step  5 (Ratio  test)  Determine 


\ . max  { max  : p±  > 0 , i€l},  max  £ -q j /p  : Pj  < 0 , J € J } }. 

If  X < 0,  terminate.  A solution  is  given  by  j • Otherwise,  let 

k " (cr,9)  be  a maximizing  arc.  Go  to  Step  1. 

According  to  the  third  remark  at  the  end  of  Section  2,  the  proposed 
algorithm  will  compute  a set  of  equilibrium  flows  in  a finite  number  of 
steps. 


Some  Computational  Remarks.  The  updating  of  the  nonbasic  components  in 
Step  4 is  easy  to  carry  out.  Indeed,  for  each  nonbasic  arc,  it  suffices 
to  determine  the  basic  arcs  which  are  adjacent  to  it.  This  provides  an 
efficient  way  to  multiply  the  product  MIJ(qJ  , pj)  from  which  the  updating 
can  be  achieved  trivially.  Furthermore,  the  relationship 


(qt  » pt)  + (qti  , Pt,)  - (ct  + ct,  , 2pt) 


34 


I • 

I 


ft 

I 


f 


for  nonbaslc  arc*  i • (ar,9)  and  i'  • (3,cr),  can  be  used  to  reduce  the 
computational  effort  In  this  step.  The  ratio  test  in  Step  5 (and  in  Step  0 
as  wall)  can  be  somewhat  simplified  by  noting  that  q^  < 0 ■>  q^,  < 0 if  i 
and  i'  are  as  just  mentioned.  This  implication  implies  that  in  carrying 
out  the  test,  at  most  one  of  the  two  arcs  i and  1*  need  to  be  considered. 
Finally,  the  various  results  established  in  this  section  guarantee  that  a 
number  of  arcs  will  not  became  basic  at  the  next  iteration.  They  can 
thus  be  Ignored  in  the  test. 
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We  have  Implemented  the  proposed  algorithm  for  solving  some  randomly 
generated  problems.  The  data  are  generated  as  follows:  The  a^  and  are 
random  numbers  lying  in  the  interval  [0,50]  and  [0,20]  respectively.  The 


cost  ctJ  is  given  by  • dt  + dj  where  d1  is  equal  to  10  if  i belongs 
to  a certain  random  index  set  and  is  a random  number  in  [0,20]  otherwise. 


reason  for  generating  the  costs  in  this  way  is  to  ensure  that  the 
triangle  inequality  will  be  satisfied.  The  code  is  written  in  FORTRAN  and 
the  runs  are  performed  on  a DEC -20  computer  at  Carnegie -Mellon  University. 
The  results  are  summarized  In  the  table  below.  The  times  are  exclusive 
of  input  and  output. 


# nodes 

# arcs 

# basic  arcs 

# pivots 

CPU  time 
(in  sec.) 

tlme/pivot 
(in  sec.) 

40 

1,760 

24 

56 

15.850 

i 1 

0.28 

80 

6,320 

45 

65 

82.456 

1.27 

120 

14,280 

70 

176 

732.421 

4.16 

160 

25,440 

108 

420 

5,682.318 

13.53 

Notice 

that  the 

number  of  arcs 

is  equal  to 

N(N-1 ) where  N 

is  the 

number  of  regions. 

Moreover,  the  number  of  basic 

arcs  can  not  exceed  N -1. 

Looking  at  the  CPU  times,  one  could  argue  that  the  proposed  algorithm 
is  perhaps  not  performing  as  efficiently  as  it  should  be.  However,  one  has 


to  realize  that  although  these  problems  could  be  considered  as  having  a 

sparse  matrix  (the  arc-arc  matrix  of  an  N-node  problem  has  2N(N-1)2  nonzero 

200 

entries,  i.e.,  density  - — l)  , they  are  still  too  big  for  the  direct 


■ 
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applies cion  of  most  currently  existing  general-purpose  quadratic  prograsming 

and/or  linear  complementarity  algorithms.  These  times  are  roughly 
-54 

proportional  to  10  xn  , a figure  which  is  really  not  big  in  comparison 
with  the  size  of  the  arc-arc  matrix.  The  performance  of  the  algorithm 
can  be  further  improved  by  taking  advantage  of  an  important  consequence 
of  the  result  concerning  the  no  transhipment  of  flows.  This  consequence 
enables  one  to  fix  at  least  half  of  the  flow  variables  at  the  value  0 and 
thus  reduces  the  size  of  the  problem  by  at  least  one  half,  it  is  stated 
and  proved  in  the  theorem  below. 

Theorem  17.  If  a^  - a^  + cag>  °»  then  the  arc  k - (a  , 0)  € J in  each 
iteration. 

Proof.  Since  there  is  no  transhipment  of  flows,  it  is  easy  to  see  chat 
Che  basis  matrix  is  nonnegative  in  each  iteration.  If  such  an  arc  k 
were  in  J at  some  iteration,  then  since  co^lementarity  is  preserved  in 
the  algorithm,  we  would  have 

which  is  a contradiction.  Q.  E.  D. 

It  should  be  pointed  out  that  this  result  has  not  been  incorporated 
in  the  coding  of  the  algorithm 

Conclusion.  In  this  paper,  we  have  proposed  and  implemented  a parametric 
linear  complementarity  approach  for  the  computation  of  equilibrium  prices 
in  an  economic  spatial  model  addressed  by  Glassey  and  others.  In  essence, 
the  proposed  algorithm  is  a specialized  version  of  the  parametric  principal 


pivoting  algorithm  applied  to  solve  the  linear  complementarity  problem 
arising  from  the  model.  Like  Glassey's  algorithm,  ours  can  be  stated 
and  implemented  based  entirely  on  the  underlying  digraph  of  the  model. 

As  mentioned  in  the  introduction,  the  model  treated  in  this  paper 
is  a simplification  of  some  other  models  of  a more  general  nature.  Zn 
the  near  future,  we  plan  to  extend  the  technique  used  here  to  treat 
various  generalizations  of  the  model  and  will  report  our  findings 
elsewhere. 
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APPENDIX 


We  shall  show  in  Che  sequel  thac  Glassey's  Algorithm  A1  [4]  fails  co 
achieve  internal  equilibrium  in  certain  situations.  Following  the  notations 
in  the  reference,  suppose  (e.,1.)  is  an  incoming  arc,  and  E • {e,,...,e  } 

and  I « {i^»...,in}  are  the  coalitions  containing  e1  and  respectively. 

Furthermore,  for  simplicity,  we  suppose  that  there  is  no  outgoing  arc. 


Thus,  for  the  new  prices  p and  p of  nodes  e,  €E  and  i,  €l 

aj  J j 


calculated  according  to  the  algorithm,  we  have 


Pe  “ pe  + 4e  * pi  “ Pi  + **  + Ce  i “ Pi 

j j j j 1 Vi  i1 


where  pfl  and  p are  respectively  the  prices  of  the  nodes  e.  €e  and  i, 61 


J J 

before  updating. 


j 


j 


Since  originally,  respective  nodes  in  E and  I have  achieved  internal 
equilibrium,  we  have  for  all  k L 


\ + and  \ + VPir>0 


which  imply 


\ + cv* ' \ + 8«  ’ - 0 "d  \ + V4  ’ \ + fli  ‘ *i  * ° 


which  in  turn  imply  that  the  new  coalition  consisting  of  E and  I will  have 
internal  equilibrium  within  themselves. 

Now  if  the  new  coalition  formed  by  EUI  attains  Internal  equilibrium, 
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Chen  we  have  for  k ■ 1 , . . . ,m  and  l » 1 , . . . n 


P.  + c j m P4  > 0 and  p.  + c,  - p >0 

*k  W h i*  e-  ~ 


Vk  % 


or  equivalently , 


P-  + ca  ± " P,  + > 0 and  p.  + c.  „ - p + 8 - 8 > 0 

*k  W * 1 *1  Vk  % 1 • “ 

But  8 “ * - Pa  + c - p . (See  [4]).  Hence 

1 e *1  Vi  *1 

(pe  + Ce.  i + P<  > _ <Pi  + Pa  + c ) > 0 

ek  Vi  *1  Ll  e1  Vi 

for  all  k - 1,...  ,m;  L ■ 1 , . . . ,n. 

It  is  now  obvious  that  the  above  inequality  does  not  necessarily  hold 
in  general,  and  thus , internal  equilibrium  is  not  guaranteed. 

Similarly,  it  can  be  shown  that  if  an  arc  goes  out  of  the  basis 
simultaneously,  then  the  same  situation  occurs. 

In  what  follows,  we  provide  a counterexample  showing  how  the  algorithm 
fails  to  maintain  internal  equilibrium. 

The  data  are  given  as  follows  (for  seven  nodes) : 

a,  - 19.481866,  a2  - 60.873993,  a3  - 10.845196,  a4  * 6.172308,  a$  - 63.791228, 

a6  * 67*071617,  ■ 37.919311;  - 7.324636  , 8,  - 0.136526,  b2  « 3.225784, 

&2  - 0.310002  , b3  - 1.884832  , &3  - 0.530554,  b4  - 7.132555  , B4  - 0.140202, 
b5  - 0.487533  , 95  - 2.051144,  bg  - 3.92217$  8g  - 0.254960,  b?  - 5.383146, 

97  ■ 0.185765  , and  ■ Ji-j|  “ cji»  *or  *U  i,j. 

We  choose  (e,i)  to  become  basic  whenever  c » min  fp,  + c - p 1 

€i  ij  V 
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Following  the  procedure  described  by  Glassey  [5],  we  have: 

Cycle  1.  Fora  coalition  ■ (4  6) 

Cycle  2.  Fora  coalition  C2  - (2  5) 

Cycle  3.  Fora  coalition  C ■ (I  2) 

Cyd*  4*  Fcra  new  coalition  C2  - (735).  Step  2c  gives  «35  - 104.497  and 

*75  ■ 2.322544.  By  algorithm,  no  arc  is  going  out. 

Cycle  5.  Fora  new  coalition  C4  - C,  UC2  - (34567).  Step  2 gives  new 

equilibrium  prices  as:  p?  - 49.516340,  pfi  - 52.516340,  p5  - 51.516340, 
P4  - 50.516340,  p3  - 49.516340.  Step  2c  gives  0?5  - 34.600461, 

035  - 109.527072  and  9^  - 8.257801  and  x45  - 2.506109.  By 
algorithm,  no  arc  is  going  out. 

Notice  that  at  this  point,  internal  equilibrium  is  not  achieved  in 
coalition  C4  - (34567) . 
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